library(MASS)

#---
#plot MISE against bandwidth for 30 subjects and 77% average probability of missing a visit
#---

load(paste("./30 Subjects/sim_pc30_250.RData", sep = ""))
pdf(paste("mise_vs_bw_30sub_", round(100 * avg.miss.prob), "prob.pdf", sep = ""), width = 5, height = 4, pointsize = 12)
matplot(x = bwSeq, y = cbind(avg.mise.exact, avg.mise.em_sl, avg.mise.em_sm, avg.mise.em_sr, avg.mise.lc, avg.mise.ll), 
	type = "l", lty = c(1,2,3,3,4,2), col = c(1,"grey",1,"grey",1,1,1), axes = FALSE, 
	xlim = c(0.25,2.25), ylim = c(0,0.8),
	xlab = "Bandwidth", ylab = "MISE", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(side = 1, at = seq(0.5, 2, by = 0.5), tick = TRUE, las = 1, cex.axis = 0.8)
axis(side = 2, at = seq(0, 0.8, by = 0.2), tick = TRUE, las = 2, cex.axis = 0.8)

min.mise <- apply(cbind(avg.mise.exact, avg.mise.em_sl, avg.mise.em_sm, avg.mise.em_sr, avg.mise.lc, avg.mise.ll), 2, min)
opt.bw <- bwSeq[apply(cbind(avg.mise.exact, avg.mise.em_sl, avg.mise.em_sm, avg.mise.em_sr, avg.mise.lc, avg.mise.ll), 2, which.min)]
dev.off()

#---
#plot MISE against probability of missing a visit for 30 subjects
#---

pdf(paste("mise_vs_prob_30sub.pdf", sep = ""), width = 5, height = 4, pointsize = 12)
tuning.par <- c("050", seq(100, 550, by = 50), seq(650, 950, by = 100))
min.mise <- matrix(NA, nrow = 4, ncol = length(tuning.par))
miss.prob <- c()
for(jj in 1:length(tuning.par))
{
	load(paste("./30 Subjects/sim_pc30_", tuning.par[jj], ".RData", sep = ""))
	min.mise[1,jj] <- min(avg.mise.exact)
	min.mise[2,jj] <- min(avg.mise.em_sm)
	min.mise[3,jj] <- min(avg.mise.lc)
	min.mise[4,jj] <- min(avg.mise.ll)

	miss.prob <- c(miss.prob, avg.miss.prob)
}
min.mise[1,] <- rep(mean(min.mise[1,]), length(tuning.par))

matplot(x = miss.prob, y = t(min.mise), type = "l", lty = c(1,3,4,2), col = 1, 
	xlim = c(0.5,0.9), ylim = c(0.05,0.2), axes = FALSE, 
	xlab = "Average Probability of Missing a Visit", ylab = "Minimum MISE", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(side = 1, at = seq(0.5, 0.9, by = 0.1), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(side = 2, at = seq(0.05, 0.2, by = 0.05), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

#---
#plot MISE against number of subjects for 77% average probability of missing a visit
#---

pdf("mise_vs_sub_77prob.pdf", width = 5, height = 4, pointsize = 12)
n.subject <- c(seq(10, 100, by = 10), 200, 500, 1000)
min.mise <- matrix(NA, nrow = 4, ncol = length(n.subject))
for(kk in 1:length(n.subject))
{
	load(paste("./", n.subject[kk], " Subjects/sim_pc", n.subject[kk], "_250.RData", sep = ""))
	min.mise[1,kk] <- min(avg.mise.exact)
	min.mise[2,kk] <- min(avg.mise.em_sm)
	min.mise[3,kk] <- min(avg.mise.lc)
	min.mise[4,kk] <- min(avg.mise.ll)
}

matplot(x = 1:length(n.subject), y = t(min.mise), type = "l", lty = c(1,3,4,2), col = 1,
	xlim = c(1,length(n.subject)), ylim = c(0,0.25), axes = FALSE,
	xlab = "Number of Subjects", ylab = "Minimum MISE", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(side = 1, at = 1:length(n.subject), labels = n.subject, tick = TRUE, las = 2, cex.axis = 0.8)
axis(side = 2, at = seq(0, 0.25, by = 0.05), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()
